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FOREWORD 

The  work  in  this  report,  was  initiated  under  Contract  DA-01-021 
AMC-11536(Z)  for  exploratory  development  of  propellants  for  missiles 
and  rockets,  and  completed  under  Contract  DAAH01-67-C-0947  for 
exploratory  development  of  solid  propulsion  technology.  Both  c on«* 
tracts  were  under  the  technical  cognizance  of  Army  Propulsion 
Laboratory  and  Center,  Research  and  Development  Directorate, 

U..  St  Army  Missile  Command. 

The  application  of  finite- element  methods  to  heat- conduction 
problems  is  an  important  way  station  to  the  successful  application  of 
these  methods  to  more  complex  time- dependent  situations— specif¬ 
ically,  to  viscoelastic  problems  of  solid  propellants  and  solid- 
propellant  rocket  motors. 

The  work  described  here  has  immediate  application  to  propel¬ 
lant  grains  and  rocket  nozzles.  But  this  method  has  general  applica¬ 
tion  beyond  solid-propulsion  technology.  Accordingly,  with  the  view 
that  broader  distribution  will  ultimately  be  authorized,  the  body  of 
the  report  contains  ho  allusion  to  propellants  pr  rockets. 
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ABSTRACT 

A  new  numerical  method  for  the  solution  of  heat  conduction 
problems  in  thermally  anisotropic,  nonhompgeneous  bodies  of  complex 
geometry  was  devised  which  is  based  on  a  discretization  concept  de¬ 
veloped  in  the  matrix  analysis  of  structures.  This  discretization 
method,  commonly  referred  to  as  the  finite  element  method,  reduces 
the  problem  formulation  to  the  solution  of  a  matrix  equation  for  the 
nodal  point  temperatures  of  the  assembly  of  finite  elements.  'The 
resulting  matrix  equation  is  stable  for  any  time  step.  The  method 
is  extremely  flexible  and  easy  to  apply.  The  metk°d  was  applied  by 
writing  a  computer  program  for  the  solution  of  heat  conduction 
problems  in  plane,  thermally  anisotropic,  nonhomogeneous  bodies. 
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Section!.  INTRODUCTION 


The  approximate  analysis  of  heat  conduction  and  other  dif¬ 
fusion  phenomena  in  bodies  of  complex  geometry  has  generally 
been  accomplished  by  using  various  finite  difference  techniques, 
e.  g. ,  [  1] .  These  methods  suffer  from  a  number  of  limitations 
or  restrictions  which  depend  on  the  type  of  formulation.  Explicit 
finite  difference  methods,  for  p sample,  have  stability  criteria  that 
often  make  the  time  increment  _  aquirements  excessively  small, 
which  in  turn  make  computation  time  excessively  large.  Regular 
grid  arrays,  v/hich  yield  simple  finite  difference  operators  are 
difficult  to  adapt  to  complex  boundaries.  This  problem  is  com¬ 
pounded  when  multi-material  bodies  are  considered,  since  each 
material  interface  must  be  treated  as  a  boundary. 

Other  types  of  solution  are  becoming  more  common*  espe¬ 
cially  those  approximate  methods  based  on  variational  principles 
[  2]  .  This  fact,  coupled  with  experience  and  ideas  developed  in 
applying  variational  methods  to  the  matrix  analysis  of  structures, 
has  led  to  the  present  development.  From  this  previous  experi¬ 
ence  it  was  expected  that  the  use  of  finite  element  methods  would 
make  multi-material  bodies  and  bodies  of  complex  geometry  more 
amenable  to  solution,  as  well  as  providing  a  compatible  nodal,  point 
system  for  coupled  usage  with  numerical  stress  analysis  proce¬ 
dures  based  on  similar  concepts. 

The  present  work  applies  a  variational  method,  along  with 
a  discretization  concept  developed  in  the  matrix  analysis  of  struc¬ 
tures,  to  numerical  analysis  of  heat  conduction  in  thermally  aniso¬ 
tropic,  nonhomogeneous  bodies.1  This  discretization  method, 
commonly  referred  to  as  the  finite  element  method,  reduces  the 
problem  formulation  to  the  solution  of  a  matrix  equation  for  the 
nodal-point  temperatures  of  the  assembly  of  finite  elements. 

I 

First,  a  functional  of  the  temperature  field  and  of  its  first 
time  derivative  is  introduced.  Then  it  is  shown  that  when  the 
functional  is  an  extremum,  it  satisfies  the  heat  conduction  equation 
throughout  the  body  and  satisfies  general  flux  boundary  conditions 
over  the  part  of  the  boundary  where  the  temperature  is  not  speci¬ 
fied.  Under  the  assumption  of  a  piecewise  linear  temperature  dis¬ 
tribution  in  a  small  quadrilateral  element  which  is  made  up  of 
four  triangular  elements  with  linear  temperature  distributions,  the 

^xter  the  initiation  of  this  work,  a  similar  approach  £0  this 
problem  was  published  by  Nickell  and  Wilson  [  3] . 


variational  principle  is  used  to  establish  a  matrix  equation  £or  the 
element’  in  terms  of  its  corner,  or  nodal-point,  temperatures  and 
its  boundary  Conditions.  Since  this  is  done  in  a  matrix  formula¬ 
tion,  the  resulting  equations  for  the  assemblage  of  finite  elements 
constituting  the ‘body  of  interest  are  easily  assembled  by  methods 
of  matrix  algebra. 

The  resulting  matrix  equation  is  stable  for  any  time  step, 
thus  offering  potential  advantages  over  the  explicit  finite  difference 
methods  in  computer  running  time.  Each  quadrilateral  element 
in  the  assemblage  may  have  distinct  and  anisotropic  thermal  pro¬ 
perties.  Complex  geometries  can  be,  approximated  as  closely  as 
desired  with  a  piecewise  linear  boundary. 

Although  the  development  is  done  in  general  terms,  the  com¬ 
puter  program  written  to  demonstrate  the  method  is  limited,  to  a 
plane,  nonhompgeneous  body  whose  axes  of  anisotropy  must  be  in 
the  same  Cartesian  frame  over  the  body.  Internal  heat  generation 
is  neglected, but  adiabatic,  constant  flux,  convective,  and  tempera¬ 
ture  boundary  conditions  may  be  applied.  Extension  of  the  program 
to  general  anisotropy,  internal  heat  generation,  and  axially  sym¬ 
metric  bodies  can  be  easily  accomplished.  Extension  to  three- 
dimensional  geometries  is  straightforward  in(  concept  but  will  in¬ 
volve  extension  of  present  programming  concepts. 


Section  H.  FORMULATION  OF  THE  VARIATIONAL  PRINCIPLE2 


Let  n,  a  functional  of  the  temperature  field  U(x,  y,  z)  and 
the  first  time  derivative  of  the  temperature  field  U(x,  y,  z),  be 
defined  by  (1). 

n(u,u) .  U^k^U.j  +  pcUtj-QuJdV 

-  f  n.-q.UdS  ,  (1) 

where 

V  =  volume  of  the  region, 

S  =  bdundarv  of  the  region,, 
k_  =  k^(x,  y,  z)  =  thermal-conductivity  tensor, 
c  5  c(x,  y,  z)  =  specific  heat, 

p  =  p(x,  y,  z)  =  density, 

Q  =  internal  heat-source  density, 

q.  =  heat  flux  vector  across  a  boundary,  and 

n.  =  unit  normal  vector. 

i 

A  comma  denotes  differentiation  with  respect  to  the  following  subt 
script,  and  repeated  subscripts  imply  summation.  The  quantities 
k,  c.  and  p  are  assumed  to  be  temperature  and  time  independent. 
Q  and  q  are  specified  functions  of  time,  and  S  and  V,  characteri¬ 
zing  the  region,  do  not  change. 

• 

The  variation  of  Il(U,  U)  with  respect  to  U  (with  U  held  con¬ 
stant)  is  given  by 


where  €  is  a  small  parameter  and  X  is  any  one  of  a  family  of 
functions  that  is  0  on  the  portions  of  S  on  which  temperature,  is 
specified  and  arbitrary  elsewhere.  An  extremum  of  the  function 
II  is  sought,  which  implies  that  6II(U,  U)  must  be  zero,  i.  e. , 

ZA  similar  variational  principle  for  isotropic  Materials 
is  given  in  [  4] . 


Starting  with 


H(U  +  ex,  U)  =  ^|i(U  +  eX),.  k„  (U  4-  eX),. 

+  pc  (U  4’  eX)TT  -  Q(U  +  €X)^  dV 

\  v  -  J  n-q^u  +  e?s.)  dS  , 


there  results 


^  =  i([(u+eX)'iVL5'kij(u+£V 

+  pcXU  -  QX^  dV  ~  n.q^  X  dS  • 


The  volume  integral 

t(U’ikijX>’j  ^ 

can  be  transformed  irf  ->  a.  surface  integral 
^(U,ikijX),j  dV  =  ^n.ki.U^XdS 
which  gives,  when  (3)  is  evaluated  at  €  =  0f 


sn(u  +  CX,  U) 
ae 


=  f  (-k..U,..  +  pcU  -  Q)  x  dV 

l/tr  3* 


e=o  "v 


+  J  n.k..U,.  Xd$  ~  J  n.q.XdS  =  0;  (4) 


The  vanishing  of  611  requires  then,  that  in  V 
kij  Uiji =  ‘ Q 


and  on  $ 


ft  8^  •  i  {&} 

Eq,  (5)  is  the  Fourior  heat  equation  and  (6)  defines  the  hast  Hus 
q  at  tho‘  surface  of  the  body.  Therefore  a  fcsctioa  U  which  gives  -fit* 
OKtremiun  of  the  functional  defined  Ly  (1)  satisfies  both  the  field 
equation  and  boundary  equations  of  transient  heat  conduction  in  ah 
anisotropic  body.’ '  ; 


tr 
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Soction  m.  £>jl£C P. E TIZATION  OF  TES  KIC3LBM 

In  the  preceding  dsvdepmcat,  in  %  funs.  .  pi  any 

functions  U  and  ft  which will  satisfy  the  boundary  conditions  on  $. 
However,  if  the  choice  of  U  and  0  is  restricted  such  that  their  coaly 
arbitrariness  is  in  certain  constants  in  their  formulation,  the  Sun?* 
tioaal  II  becomes  a  real-valued  function.  In  particular,  it tho  ccs&. 
stants  are  the  vector  of  nodal-point  values1,  jj  and  6  of  U 
and  ft,  II  becomes  H(}j,  j$).  Finding  an  extremum  of -this  real- 
valued  function  is  equivalent  to  satisfying  thp  following. 

8n{a,$ 


In  the  following,  the  body  will  bo  considered  to  bo  divided  in¬ 
to  a  number  of  tetrahedral  or  piano  triangular  elements.  Those  are, 
in  some  sense,  small  with  respect  to  the  temperature  gradient  and 
boundary  contours  such  that  the  temperature  distribution  and  boundary 
can  be  represented  by  a  piecewise  linear  approximation.  The  nodal 
points  for  the  numerical  analysis  will  be  the  vertices  of  tho  elements. 


Let  tke  temperature  field  in  an  element  be  given  by 

U(x,  y,z,t)  =  $ (x,  y,  a)  £»(t)  (3) 

and  the  time  rate  of  temperature  change  in  the  element  be  given  by 

U(x,  y,  z,t)  *  $x,  y,z)Bu  (t)  ,  (9) 

where  <f>  and  >Jj  are  vectors  which  specify  the  spatial  dependence  of  U 
and  ft  and  and  &  are  the  vectors  of  nodal  point  values.4  The  matrices 
of  constants,  ^  and  ||,  are  defined  by  the  above  relationships. 

The  temperature  gradient  U,,  in  the  element  can  be  expressed 
in  terms  of  the  nodal-point  temperature^  by 


=  Au  •  (10) 


AU  = 


U, 

u, 

U, 


x 


2 


3The  notation  4  indicates  the  matrix  A,  and  jj  the  vector  u. 

^Xt  should  be  noted  that  throughout  this  development,  the  fields 
U  and  ft  have  been  taken  to  be  independent.  In  the  computer  program, 
however,  $  andjjj,  and  therefore  A  and  §,  were  taken tobe the  same. 
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Y/riting  Q  in  terms  of  nodal  point  Quantities, 

*PC&  £  g  &  S'  Q/^ 


r  t  at  ,  t 
•  \  S  *1  23, as  . 

(ii) 

h-f) 

■ 

Taking  the  first  variation  with  respect  to  j 
it  equal  to  zero,  there  results 

and  setting 

6n  =  Ku  »  Q  #  +  C  u  -  q#  s  0  , 

(12) 

where  P  T  T 

K  -  \  A  <b*  kd>*  AdV 

a  Jy  S  &  B  S3  B 

(13) 

S  -  £  Pc  aVt  1  ■  • 

(14) 

9*  =  £ 

V 

(15) 

and  .  _  _ 

q*  *  f  A  <t>  Sids 

(15) 

Boundary  conditions  of  four  types  will  now  be  considered; 

(1)  Specified  temperature,  u,  =  constant  (boundary  segment 

Si), 

(2)  Specified  flux,  q  =  q  (boundary  segment  S2), 

(3)  Convective,  q  =  h(u.-u<j),  where  h  is  a  film  coefficient 
and  uq  is  the  environmental  temperature  (boundary  seg¬ 
ment  S3),  and 

(4)  Adiabatic,  q  =  0  (boundary  segment 

The  boundary  integral  (16)  now  becomes 

Jl*  =  q*+Hu-h*  ,  (1?) 

*•*  ** 


7 


whore 


H 

53 


=  h 


* 

Ss 


(18) 

0?) 


(20) 


The  integral  over  Sx  is  zero  since  the  variation  of  the  functional  was 
specified  as  zero  over  that  portion  of  the  boundary,  and  the?  integral 
over  $4  is  zero  since  there  is  no  heat  flow  across  the  boundary. 

To  assure  the  extremum  of  the  functional  II,  it  is  necessary 
then  to  find  the  nodal-point  temperatures  which  satisfy  the  following 
matrix  equation. 

(K  -  H)  u  +  C  u  -  q*  -  Q  *+  h*=  0  .  (21) 


V 


Section  IV.  SOLUTION  OF  THE  GOVERNING  MATRIX  EQUATION 
I.  Solution  Method, 

To  solve  (21),  note  that  &  and  ji  are  functions  o£  time  and  Q&, 
g*  and  h#  are  given  functions  of  time.  Let  the  time  variable  bo  xe*» 
strietsd  to  the  following  set  of  variables. 


t.  =  i  A  t  ,  i  =  0, 1,  2,  „ . . 


Subscripts  ni"  in  the  subsequent  development  indicate  that 
the  subscripted  quantities  are  evaluated  at  t  =  t^. 


Let  (21)  be  written  as 
K  u.  -S-  C  u.  =  f.  , 

a  ~i  »  ~i 


where 

and 


K  =  K  -  H 

a  s  » 

f.  =  Q*  +  q*  -  h*  . 

f**X  *+  *3  m 


If  u  is  assumed  constant  for  t^  £  t  S  t  then  5^  4 
[  -  u.]/  At  and  from  Taylor*  s  expansion  about  t=  j 


(22) 


=  u. 

+ 

At 

• 

+ 

At2 

#• 

u. 

i 

2 

~  i 

=  u. 

+ 

At 

• 

u. 

+ 

At 

C  “ 

2 

A#] 

=  u. 
~1 

+ 

At 

2 

[ 

“•0.1  +  “•! 
<ri+l  ~l 

*  .  -1 

i+ 1  ~r 


(23) 


Eqs.  (22)  and  (23)  now  are  sufficient  to  determine  u.  .  and  u.  , 
in  terms  of  u^andu..  Solving  (23)  for  u^j  yields  1  "" 5 


u. .  , 

~Z+ 1 


2 

At 


[Vi 


-  u.l 


u.  . 
~1 


(24) 


Substituting  this  value  into  (22), gives 


$  *  aT 


cl 


;)ui+l  =  ^4*  +  u.l 

/  ~i+ 1  ~i+l  »L,At  ^i] 


(25) 


Also  froth  (22), 


v 


Substituting  (26)  into  (25)  results  in 

At  +  -fr  C }  u. , .  =  (■—  C  -  k\  u.  +1  *  . 

At  fsj  *44*1  \  At  a  a  J  *•  i  «-i  -4+1 


A  simpler  computation  results  by  rewriting  (27)  as 


where  u.^j  is  found  from  the  auxiliary  calculation® 


(26) 

(37) 

(38) 

(29) 

(30) 


For  the  solution  of  the  heat  flow  in  a  multi-element  body,  it  is  nec¬ 
essary  to  assemble  the  element  matrix  equations  (29)  into  a  single 
matrix  equation.  This  assembly  is  a  complex  task  which  can  be 
performed  in  an  efficient  manner  by  a  computer.  The  general  method 
of  assembly  for  matrix  equations  is  given  in  [5] ,  Section  7.  2. 

2.  Stability  of  Solution  Technique 


To  study  the  stability  of  the  solution  technique  defined  above 
by  (27),  i.e. ,  the  effect  on  the  numerical  solution  of  an  error  intro¬ 
duced  at  some  step,  consider  the  vector  jj.  which  satisfies  exactly 
the  relation 


_2_ 

At 


u. 

-I 


=  f.  +  f.AI 

~i  »»i*f  1 


(31) 


I* 


J 


i 


5Also  note  that  if  At  is  very  large  (29)  reduces  to  K  U.  =  f  *  ■ 
the  steady-state  Torm.  Thus  the  steady-state  solution  catfbeVbtajned 
in  one  iteration  from  the  computer  program  simply  by  making  the  sin¬ 
gle  time  increment  very  large. 
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Suppose  that,  at  some  stop  (i+1)  in  the  calculation  of  jj,  an  error  . 

(say  round-off  error)  is  introduced  in  the  calculation  (31),  which 
can  be  rewritten  in  terms  of  the  incorrect  value  u*  as 

<?+lt  §>h?+i  +  <§-Jt  §>%*  =  h+k*  •  <32> 

Then  subtracting  (31)  from  (32)  results  in  a  recursive  relation  (33) 
for  the  error  in  u  at  step  N.  (N  >  i+1).  Letovs  u^  -  ;  then 

&  +  g>SN  +  (S-|t  S>Vl=2  •  '  N  »  M2 .  <33> 

Solving  (33)  for  results  in 

e.T  =  CK  +  TT  G)"1  (-K  +  C)  e.T  .  =  A  e„  .  . 

~N  w  At  a'  '  At  »'  ~N-1  Jsr^N-1 

It  follows  inductively  that 

%  =  A  >  n  =  N-(i+l) 

where  e  is  the  error  introduced'  at  n  =  0.  Let  X*  be  the  absolute  value 
of  the  largest  element  of  the  mXm  matrix  A,.  Then 

(mX*  e  >  A  e  . 

Consider  solutions  of  the  form 

®N  =  (mX*)Ne  =  XH  e  ,  n  =  N  -  (i+l)  (34) 

where  X  is  a  positive  constant.  The  error  so  defined  is  greater 
than  or  equal  to  the  true  error.  Substituting  (34)  into  (33), 

§I  +  S  §>  »ls3  g[|  t-2  > 

or,  rearranging  things  slightly,  ’ 

c  (1  -  xfj  e=0_  .  •  (35) 

Defining 


11 


(35)  can  be  written 


K  -m 


(36) 


If  C  is  a  positive  definite  matrix,  then  according  to  Wilkinson  [  6,  p. 
34P  the  eigenvalues  <o  are  all  positive*  This,  in  turn,  requires  that 
|  X|  <  1*  It  follows  from  (34)  that  the  error  will  decrease  as  N  in¬ 
creases.  Since  |  X|  <  1  for  anjr  value  of  At  >  0,  the  solution  scheme 
is  unconditionally  stable  if  C  is  positive  definite.  This  property  of  C 
is  dependent  ontheassumedforms  for which,  as  indicated  in  Section 
V.  1,  give  a  positive  definite  C  for  the  forms  assumed  in  the  present 
development. 
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♦ 


the  matrix  from  its  definition  in  (37)  and  (39),  can  be  £ov  d  to  be 

0 


A  = 

"  ajbk"akbj 


ajV\bj  0 


b.-b. 

1  k 

L  Vaj 


bk  -bj 


-a,. 


a. 

J 


(41) 


BqSi  (37) -(41)  are  sufficient  to  define  a  linear,  spatial,  temperature 
field  and  a  linear,  spatial,  temperature-rate  field  in  terms  of  the 
nodal  point  values  of  the  temperature  and  temperature  rate,  respec¬ 
tively. 


Since  A  and  dM  are  constant  over  the  element,  (13)-(16)  may  be 
written  a 


if  K,  pc;  and  Q  are  also  taken  as  constant  in  an  element. 

RS 

These  integrals  are  easily  evaluated  in  terms  of  the  first  and 
second  moments  of  area. 

The  boundary  integrals  in  (l8)-(20)  also  simplify  to 
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The  coefficient  matrices  appearing  in  (29)  may  now  b©  written 


K-  ATd>»T  k<fc>A  A  -  h  A 

PS  (3  «  fit  K  »  « 


■'IX 


d$  A 


c.pc/py*#]  a  , 


i  =  QAT  r^TdA/>{-  ATJ  ^Tn  J  ds  -  hu0AT  f  <pT4 $ 

=  AT  £q  J  $T  dA  4  J  ^Tn  3  dS  -  hu0  £  .(51) 

A  S  2  §3 

In  the  development  of  the  computer  program  <p  was  taken  in  the  linear 
form 

<£(*.y)  =  {i.x.y}  , 

In  this  case,  for  a  triangular  element,  (50)  becomes 


2  1  1 
1  2  1 
1  1  2 


which  is  positive  definite  as  required  for  stability  in  Section  IV,.  2. 

2.  Quadrilateral  Element  Matrices 

It  is  convenient  in  terms  of  programming  logic  to  work  with  a 
quadrilateral  element.  For  this  purpose  a  quadrilateral  element  com¬ 
posed  of  four  triangular  elements,  as  shown  in  Fig.  2,  was  used  in 
the  present  computer  program.  The  four  triangles  are  determined 
by  defining  the  coordinates  of  the  common  point  to  be  the  average  of 
the  coordinates  of  the  other  four  points.  The  common  point  is  elimi¬ 
nated  from  explicit  representation  by  the  following  procedure. 
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ETG-  JL'.  .  QUADRILATERAL.  ELEMENT 

The  matrix  equation  for  this  quadrilateral  element  can  be 
expressed  by  assembling  the  matrix  equations  for  each  triangular 
element  by  addition  of  terms  at  each  nodal  point  in  the  manner 
used  in  the  direct  stiffness  method  of  structural  analysis.6  Eq. 
(22)  for  a  quadrilateral  element  takes  the  form 


Ka  KJ2  K23  K15 

Xji  K22  Kjj  K 24  K05 

K3J  K32  K33  K34  K35 

^41  K*2  K43  Kji4  K45 

KS1  Ksa  Ks3  Kst  Kgs 


’ll  32  '-'13 


0  C« 


u2.  C2I  C22  0  C24  C25 

u3  +  C31  0  C33  C34  C35 

U4  0  C42  C43  C44  G45 

u3  Cgj  C52  C53  G54  C55 


6See  [5] ,  Section  7.2  for  a  description  of  this  assembly  method 


Bewriting  this,  there  results 


4 


m 


[Kj^l  [KjS] 

‘[Ujf 

4- 

lcJ 

r*  -oi 

inf 

* 

& 

jK5j]  [Kb]_ 

Gss- 

- 

,l,j=l,2,3,4  , 
(53) 


4.  »  •  v*  xpv  ■»- 

tors.  The  subscripts  i,>  j  now  represent  .nodal  points,  instead  o£ 
time  increments.  Eq.  (53)  can  then  be  written  as  two  equations, 


and 


+  [K^3u^+  [Ci.][u.]  4-  [CjJus 


■ lt? 


(54) 


[Kg.][u.]  4-  .K55US  4*  £,Csi]:[uJ  4-  Cg5  ug  —  is, 
J.  .  J  J  ) 


(55) 


The  interior  nodal  point  quantities  ug  and  65  cannot  be  elimi¬ 
nated  from  (54)  by  use  of  (55)  as  it  stands.  However,  if  thfe  specific 
heat  matrix  C(5  x  5)  is  approximated  by  lumping  the  heat  capacitive 
at  the  four  external  nodal  points,  C  becomes  a  diagonal  matrix*  with 
G53  =  0  and  (54)  can  be  written 


[K.jJLuj]  4-  [K.s]us  4-  =  [f.].  (56) 

and 

[K5jj  [  Uj]  4-  Ksg  Ug  *  fg  ,  (57) 


in  which  [G. .]  is  now  the  (4  x  4)  submatrix  of  the  diagonal-lumped 
specific  heat?  matrix. 

Solving  (57)  for  ug  and  substituting  into  (56),  there  results 
[Kj.Hu.]  +  [K.5]K55‘l{fs  .  [K5j][u.]}  +  [Cj.Hu.]  =[£j  (58) 
or 

{[Kj.]  -[K15]KI,'1[K3j]}  [u.J  +  [G.jl [Uj]  =[£.]HKi5lKss'ii5.(59) 


✓ 


1 


This  form  of  the  specific -heat  matrix  is  also  positive 


definite. 
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This  equation,  (59),  au.~  is  analogous  to  (22)  except  that  Etand.C  aye 
now  4X4  matrices  and,  £  is  a  4X  i  vector,.  , . 


1 


Section  VI.  COMPUTER  PROGRAM 
Description 

The  organization  and  coding  of  the  present  computer  pro¬ 
grams  rely  heavily  on  concepts  developed  previously  in  finite 
element  structural  analysis  programs,  particularly  those  described 
in  f  5]  .  Two  programs  are  described  below.  AMGG42  is  the  heat- 
conduction  program,  and  AMG042P  i 3  an  associat'ed  plot  program 
which  may  be  used  to  aid  in  reducing  the  output  data  to  graphical 
form. 


Program  AMG042,  which  has  been  written  to  effect  the  solu¬ 
tion  of  the  matrix  equations  formulated  in  Section  V,  is  somewhat 
more  restricted  than  that  development.  Although  the  steps  necessary 
to  generalize  the  program  are  obvious,  these  are  not  necessarily 
trivial.  Presently  the  directions  of anisotropy  of  conductivity  of  each 
element  must  all  lie  in  the  same  Cartesian  frame.  Likewise  there 
is  no  provision  for  internal  heat  generation.  However  the  material 
properties  may  vary  from  element  to  element. 

The  development  has,  in  general,  been  applicable  to  bodies 
of  fairly  arbitrary  shape.  However,  the  necessity  of  employing  a.  for¬ 
mal  solution  method  consistent  with  minimum  effort  in  data  input 
has  resulted  in  some  restraints  in  the  computer  program.-.  The  net¬ 
work  of  quadrilaterals  needed  for  solution  was  regularized,  with  a 
two-dimensional  nodal-point  identification  array.,  which  then  pro¬ 
vided  a  systematic  framework  for  solution  of  the  matrix  equation. 

This  grid  method  was  first  developed  for  stress-  analysis  purposes, 
and,  although  it  is  described  in  some  detail  in  Section  VI.  2,  a  more 
comprehensive  treatment  is  given  in  [  6],  Aside  from  the  require¬ 
ments  on  grid  network,  some  further  restrictions  are  imposed  by 
the  boundary  condition  subroutines  which  are  described  below. 

In.  setting  up  the  program  logic,  it  became  obvious  that  in¬ 
cluding  completely  general  time-dependent  boundary-condition  op¬ 
tions  for  arbitrary  geometry  would  be  extremely  difficult.  There¬ 
fore,  it  was  decided  to  handle  the  boundary  conditions  by  separate 
short  routines  to  be  prepared  for  each  class  of  problems.  The 
boundary  condition  subroutines  included  in  this  report  are  written 
to  apply  only  to  a  rectangular  nodal-point  identification  array.  This 
does  not  imply  that  the  program  in  its  present  form  is  limited  to  a 
rectangular  region. 
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The  sequence  of  operations  of  AMG042  is  given  by  the  flow 
chart  shown  in  Fig.  3.  The  coefficients  of  the  complete  matrix 
equation  are  assembled  from  the  coefficients  of  each  quadrilateral 
in  a  manner  analogous  to  the  direct  stiffness  method  of  structural 
analysis.  See  [  6]  t.  p.  28,  for  a  more  detailed  description  of  the 
assembly  process.  Modifications  for  boundary  conditions  are  made 
in  a  similar  manner. 

2.  Mesh  Layout  and  Generation 

l 

The  requirement  of  closely  approximating  the  contours  of 
complex  regions,  together  with  the  desirability  of  a  fine  mesh  size 
and  its  attendant  high  accuracy,  makes  the  use  of  a.  large  number  of 
nodal  points  desirable.  The  program  allows  the  user  to  employ  a 
maximum  of  496  quadrilateral  nodal  points.  Obviously  the  layout  and 
specification  for  the  program  of  the  locations  of  such  a  number  of 
points  is  a  tedious  and  time-consuming  job  in  which  the  probability 
of  human  error  is  high.  To  minimize  this  effort  and  to  preserve  as 
much  general  utility  as  possible,  a  scheme  for  the  internal  (to  the 
program)  generation  of  much  of  the  required  data  has  been  incorpo¬ 
rated  in  the  program*  This  same  scheme  has  been  used  previously 
in  stress  analysis  programs  [  6]  Certain  restrictions  are  imposed 
on  the  layout  of  the  nodal  points ,  but  the  reduction  in  the  effort  re¬ 
quired  to  effect  the  solution  of  a  given  problem  adequately  compen¬ 
sates  for  these  restrictions. 


To  lay  out  a  nodal-point  system  for  the  body  to  be  analyzed, 
the  region  of  the  x-y  plane  constituting  the  body  is  covered  (insofar 
as  any  curved  boundaries  will  permit)  with  an  array  of  convex  quadri¬ 
laterals.®  Each  vertex  of  a  quadrilateral  is  called  a  nodal  point  or 
node.  Each  nodal  point  is  identified  by  an  ordered  pair  of  positive 
integers,  denoted  by  (l,  j).  The  nodes  may  thus  be  thought  of  as  a 
subset  of  the  lattice  points  in  the  I-J  plane.  Nodes  with  common 
second  member  J  are  said  to  lie  in  the  same  row,  although  this 
implies  nothing  about  their  location  in  the  x-y  plane. 


The  scheme  for  mesh  generation  may  be  thought  of  as 
representing  a  one-to-one  mapping  from  the  I-J  plane  into  the  x-y 
plane.  Fig.  4  illustrates  this  mapping.  The  points  in  the  I-J  plane 

^Tbe  use  of  a  quadrilateral  element  with  a  vertex  angle 
greater  than  180°  may  result  in  erroneous  calculations  for  that 
element.  A  vertex  angle  of  180°,  which  is  acceptable,  gives  the 
quadrilateral  the  appearance  of  a  triangle. 
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FIG.  4.  DEMONSTRATION  GRID 


i  i  ; 


$ 

% 


22 


are  shown  in  Fig.  4a  and  their  image  points  in  the  x-y  plane  are  shown 
in  Fig.  4b.  It  can  be  seen  that  the  inverse,  images  of  the  quadrilat¬ 
erals  in  the  x-y  plane  are.  squares  in  the  I-J  plane*  Each  quadrilateral 
(or  continuum  element)is  identified  by  the  I.,  I  coordinates  of  the  node 
whose  inverse  image  lies  at  the  lower  left-hand  vertex  of  the  inverse 
image  of  tne  quadrilateral.  Thus  the  nodes  which  are  vertices  of 
element  (I,  J)  are  the  nodes  (I.  ,T).  (l+l,j)r(l,j+l),  and  (Hi,  J+T). 

It  may  be  noted  that  not  every  boundary  node  need  have  an  element  , 
associated  with  it.  In  Eig.  4*  circles  represent  nodes  associated 
with  elements  and  squares  those  which  are  not  associated  with  any 
element.  The  unfilled  circles  represent  nodes;  whose  coordinates 
Were  generated  by  the  program. 

I 

An  important  restriction,  which  is  due  to  the  bookkeeping 
procedure  used  in  the  program  toasembie  the  element  stiffness  into 
the  stiffness  for  the  entire  structure,  may  be  phrased  thus:  if,  in 
any  given  row,  IMIN  and  IMAX  are  respectively  the  least  and  greatest 
value  of  I  for  which  there  is  a  node,  then  there  must  be  a  node  in 
that  row  for  each  I  such  that  IMIN  2I<  IMAX.  For  the  present 
program  IMAX  :£  16  and  JMAX  «S  31,  The  limiting  values  of  IMAX 
and  JMAX  may  be  va_ied  by  changing  the  appropriate  dimensions  in 
the  COMMON  statements  so  as  to  stay  within  the  capacity  of  the  com¬ 
puter.  All  nodal  points  that  define  the  boundary  must  have  their 
coordinates  specified, and  any  other  nodal  points  may  either  be 
specified  or  calculated  by  the  internal  generation  scheme. 

The  mesh  generation  is  accomplished  in  the  following  manner. 

A  data  card  containing  the  values  of  I,  J  and  the  x,  y  coordinates  is 
input  to  the  computer  for  each  node  Whose  coordinates  are  to  be 
specified.  Such  nodes  must  include  at  least  all  nodes  on  the  bound¬ 
ary  of  the  region  of  interest,  as  well  as  on  any  interfaces  between 
regions  of  different  materials.  As  many  other  nodes  as  the  u3er 
may  desire  may  have  their  coordinates  specified,  but  no  others  are 
necessary.  As  the  data  cards  are  read,  a  list  is  cor-.piled  of  the 
minimum  and  maximum  values  of  I  for  each  J,  and  each  node  for 
which  coordinates  have  been  input  is  identified  and  the  coordinates 
are  stored. 

An  option  is  included  to  permit  the  input  of  straight-line 
segments,  corresponding  to. I  =  constant  or  J  =  constant,  which 
are  to  be  divided  into  equal  increments.  The  I,  J  corresponding  to 
the  smallest  I  (or  smallest  j)  is  input  in  the  first  position  on  the 
card,  with  the  I,  J  corresponding  to  the  largest  I  (or  largest  J)  being 
input  in  the  secopd  position.  Corresponding  x,  y  coordinates  ate  in¬ 
put  into  the  first  and  second  coordinate  positions.  The  line  segment 


is  inter  nally  divide  u.  and  assigned  equally  spaced  nodal  points*  Note 
that  TYPE,  BCCOBE,  and  IJCGDE  (doseribsdin  AppOndfeaA, .  l.a) 
must  be  the  same  tor  ail-nodal  points*  If  only  a  single  code  is  to  Be 
input,  the  second  I  and  J  positions  and  the  second  coordinate  positions 
are  left  blank.  A  polar-coordinate  input  option  is  also  provided. 

After  all  the  desired  nodal  point  cards  have  been  input,  the 
coordinates  for  all  unspecified  nodes  which  have  I  in  the  interval 
IMIN  <  l  <  IMAX  for  the  proper  J,  are  calculated  for  all  J,  The 
calculation,  or  mapping,  of  the  coordinates  ip  achieved  by  solving 
twice  the  finite-difference  analogue  of  Laplace1  s  equation  on  the 
lattice  points  in  the  I-J  plane.  First,  the  x  coordinates  of  the  bound¬ 
ary  points  are  used  as  boundary  values  of  the  unknown  harmonic 
function,  and  the  functional  values  obtained  on  the  interior  points  are 
taken  as  the  x-coordinates  of  the  corresponding  injage  points  in  the 
x-y  plane*  A  similar  procedure  yields  the  y-epordinates  of  the  un¬ 
specified  nodes.  It  should  be  noted  that,  in  general,  this  method  tends 
to  yield  nodal  points  with  uniform  spacing.  If  this  is  not  deemed 
desirable,  some  nodal  points  interior  to  the  region  may  have  their 
coordinates  specified  to  control  the  distribution  of  the  remaining 
points.  ' 


Section  VII.-  .ILLUSTRATIVE  PROBLEMS 

Several  problems  to  illustrate  the  utility  and  accuracy  of 
the  program  have  been  solved  and,  when  possible,  compared  with 
formal  solutions.  Since  these  are  of  the  form  of  illustrations,,  the 
problems  are  posed  in  dimensionless  form  wherein  any  consistent 
set  of  units  may  be  inferred.  Unless  otherwise  stated,  K,  p,  and 
c  were  taken  to  be  unity.  *  ! 


To  demonstrate  the  accuracy  of  the  solution  technique,  th^ 
problem  of  an  isotropic,  homogeneous  2X  2  square  initially  at  a  uni¬ 
form  temperature  of  1  with  boundaries  held  at  0  was  solved.  For 
a  one-quarter  symmetric  section  of  the  square,  a  14 X  14  grid  was 
used.  The  resulting  temperature  distributions  along  three  constant 
coordinate  lines  are  shown  in  Fig.  5  and  compared  with  theoretical 
results  from  Carslaw  and  Jaeger  [  7]  .  Agreement  is  quite  good. 


C  Onvective,  Boundary  C  ondition 


The  problem  of  a  hollow,  circular  cylinder  with  convective 
boundary  conditions  was  run.  Because  of  the  assumed  symmetry 
it  was  only  necessary  to  run  a  sector-shape  geometry  with  adiabatic 
boundary  conditions  on  the  straight  sides  and  convective  conditions 
on  the  inner  and  outer  boundaries.  A  45®  sector  was  used,  with  an 
inner  radius  of  0.25  and  outer  radius  of  1.  The  convective  boundary 
conditions  h  (T-T^  used  were  35.0  on  the  inner  boundary  and  70.  0 
on  the  outer  boundary.  Initial  temperature  was  zero,  and  the  enviro¬ 
nment  temperature  was  1. 


Fig.  6  illustrates  the  comparison  with  the  results  from  a 
finite  ^difference  program  [$] ,  As  can  be.  seen,  agreement  is 
essentially  perfect.  For  this  particular  run,  the  time  increment 
for  the  finite  element  solution  was  taken  as  0.001  while  the  time 
increment  was  0.000125  for  the  finite  difference  solution.  When 
the  time  increment  for  the  finite  element  solution  was  taken  10  times 
larger,  0.01,  the  oscillations  shown  in  Fig.  7  occurred.  Note  that 
the  boundary  temperature,  as  indicated  frojn  results  with  smaller 
time  increments,  should  have  reached  over  90%  of  its  total  tempera¬ 
ture  change  during  the  first  time  increment  of  0.01,  Despite  this 
crudeness  and  the  resulting  oscillations  near  the  boundary,  the 
solution  near  the  center  of  the  slab  is  fairly  accurate  for  all  times 
and  the  solution  near  the  boundaries  becomes  more  accurate  for 
increasing  time  as  the  oscillation  dies  away.  This  is  illustrated 
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FIG.  5.  TEMPERATURE)  IN  A  SQUARE 


in  Fig.  8  in  which  the  temperature  at  r  =  0. 25,  r  =  0.60,  and 
r  =  1.0  a-ce  plotted  versus  time. 


3.  Flux.  Boundary  Condition 

The  behavior  of  constant-' flux  boundary  conditions  were  in¬ 
vestigated  for  a  rectangular  slab  witha  constant  flux  and  su’ppliedafrfr 
opposite  faces  while  the  other  two  faces  were  adiabatic.  The  slab 
was  initially  at  a  temperature  of  zero*  The  results  are  illustrated 
in  Fig.  9.  Only  one-half  of  the  slab  is  illustrated*  The  center  is 
on  the  left  of  the  figure.  The  lines  are  from  the  series1  solution  of 
Carsla,w  and  Jaeger  [•?] .  * 

4 0  Nonhomogeneous  Properties' 

An  axisymmetric  cylinder  with  conductivity  and  specific 
heat  which  vary  inversely  with  radius  was  studied,  Ipitial  tempera¬ 
ture  of  the  cylinder  was  given  as  zer«%  and  the  internal  and  external 
boundary  were  subjected  to  a  temperature  of  1  at  time  t  =  0.  The 
results  are  compared  in  Fig.  10  with  the  formal  solution. 

If  K  =  —j?—  and  pc  =  ,  then  for  an  axisymmetric  cylinder 

the  heat-conduction  equation  becomes 


which  is  the  equation  for  a  homogeneous  slab.  This  solution  was 
obtained  from  [  7,  p,  10 1]  to  plot  in  Fig.  10.  Kq  was  taken  to  be 
unity  and  the  product  p0o  was  taken  to  be  5.  Agreement  with  the 
formal  solution  is  very  good  despite  the  crude  mesh  of  nine  radial 
increments. 


5,  Anis otr opic  C  onductivity 

The  quadrilateral  shown  in  Fig.  II,  with  the  conductivity  in 
the  x-direction  equal  to  4  times  the  conductivity  in  the  y-direction, 
was  used  to  check  the  anisotropic  features  of  the  program.  This 
can  be  checked  with  an  isotropic  body  by  the  following  analogy. 

Let  kx,  k  ,  x,  and  y  represent  the  conductivity,  and  coordinates 
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FIG.  9.  TEMPERATURE  IN  A  SLAB  Y/ITH  CONSTANT-FLUX 
HEAT  INPUT  AT  THE  BOUNDARY 
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If  we  define  a  new  coordinate  s£  =  x,  the  equation  may  be  written  as 


82u 


8su 


Le^2  '  ay2  J 


pc 


du 

dt 


Thus,  the  solution  of  an  isotropic  problem  in  x,  y  coordinates  with 
isotropic  conductivity  k„  gives  a  temperature  field  similar  to  an 
anisotropic  problem  with  =  n2k^  and  x  =  nx. 

For  the  present  problem,  the  quadrilateral  is,  in  fact,  a 
square  with  the  x  coordinate  doubled,  and  with  =  4  ky.  For  initial 
conditions  of  T  =  0  with  the  boundaries  held  atT  =  1,  the  transient- 
conduction  problem  was  worked  for  both. the  isotropic  square  and  the 
anisotropic  quadrilateral.  The  temperature  calculated  in  the  two 
problems  agreed  to  five  significant  figures.  Temperature  contours 
for  the  quadrilateral  are  shown  in  Fig,  H  for  time  t  =  0,  05,  =  4, 

ky  =  1.0,  pc  =  0, 16,  A  mesh  of  16  x  31  was  used, 

6.  Complex  Geometry 

As  an  example  of  the  utility  of  the  program,  an  example  is 
given  in  Figs.  12  through  15  which  demonstrates  its  use  on  a  one- 
sixth  symmetric  section  of  a  cylinder  with  a  star-shaped  perfora¬ 
tion  subjected  to  severe  convective  cooling  conditions.  The  geometry 
with  the  internally  generated  finite  element  grid  is  . shown  in  Fig, 

12.  The  initial  temperature  of  the  body  was  T^  and  the  environment 
temperature  TQ.  Results  are  presented  in  c'  nensionlsss  quantities. 
Isotherms  are  demonstrated  in  Fig.  13,  and  In  Figs.  14  and  15, 
temperature  profiles  are  compared  with  those  of  Willoughby  [9j. 
Willoughby' s  solution  is  shown  in  solid  lines.  Willoughby  used  a 
combination  of  conformal  mapping  and  finite  differences  to  obtain 
his  solution. 
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12.  GRID  LAYOUT  ON  SYMMETRIC  SECTION  OF  SIX- POINT  ST  AR~  PERFORATED 
CYLINDER 


V-i 


i  > 


f  ' 
>  /' 


T  *  TEMPERATURE  t  «TIM£ 

Ti  =  INITIAL  R  *  OUTER  RADIUS 

TEMPERATURE  t  .  j£JL 

Te- ENVIRONMENT  r2 

TEMPERATURE 

-  THERMAL 
D1FRJSIVITY 


| , 
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t  =o 


FIG.  14.  THERMAL  HISTORY  COMPARISON  AT  SECTION  A- A 
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3? 


TEMPERATURE 


FIG.  15.  THERMAL  HISTORY  COMPARISON  AT  SECTION  B.B 
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APPENDIX 


PROGRAM  OPERATING  INSTRUCTIONS 
1,  Finite  Element  Heat  Conduction  Program  AMGQ42 
a.  input  Data 

Initial  temperatures  mast  be  specified  at  every 
point  by  an  appropriately  coded  subroutine.  The  present  subroutine 
is  only  usable  for  a  uniform  initial  temperature  in  the  body. 

Constant  temperature,  constant  flux,  convection, 
or  adiabatic  conditions  may  be  specified  at  any  of  the  boundary 
points.  Subroutines  BCTEMP  (constant  temperature),  BCCOND 
(constant  flux),  BCCONV  (convection)  for  applying  a  boundary  condi¬ 
tion  on  any  of  four  sides  *.£  a  rectangular  nodal  point  array  are  listed 
in  this  report.  The  adiabatic  condition  is  imposed  by  the  absence 
of  other  boundary  conditions. 

In  the  present  program,  as  many  as  twenty  sets  of 
material  properties  may  be  specified  and  assigned  to  arbitrary 
elements. 


An  element  is  identified  by  the  smallest!  and  smallest 
J  associated  with  the  nodal  points  which  are  its  vertices  and  is  said 
to  be  associated  with  the  node  which  is  also  identified  by  this  I  and  J, 
The  kind  of  element  associated  with  a  nodal  point  is  specified  by  a 
three-digit  symbolic  word  TYPE  (see  Table  I).  The.  last  two  digits 
specify  the  set  of  thermal  properties  for  the  element.  The  first 
digit  indicates  whether  or  not  the  coordinates  of  the  nodal  point  are 
specified  and  if  an  element  is  associated  with  that  nodal  point. 

Two  additional  codes  are  used  for  each  nodal  point 
to  g  pertinent  information  relating  to  each  element,  BCCODE 
(see  ...  J„e  II)  branches  the  program  to  the  correct  boundary- 
condition  subroutine.  IJCODE  (sea  Table  III)  indicates  the  nodal- 
point  line  oegmeat  to  which  the  boundary  condition  applies.  Only' 
one  segment  per  nodal  poi  may  be  specified.  The  words  TYPE, 
BCCODE,  and  IJCODE  are  combined  internally  into  a  single  word, 
CODE,  to  conserve  storage  locations,  *  Of-"  is  output,  for  checking 
purposes,  with  the  coordinates  of  that  node.  The  first  three  digits 
of  CODE  are  TYPE,  the  fourth  digit  is  BCCODE  and  the  fifth  digit  i3 
IJCODE'. 


Table  I  Value  and  Meaning  of  Symbolic  Word  TYPE 

_ _ _  imi  ii ■  ~i 1 1 Miir— nrr i — — 

Value 

Meaning 

X01,  JX02,...*,X26 

Identifies  the  particular  set  of  material 
properties  to  be  associated  with  the  ele¬ 
ment.  Type  X00  is  equivalent  to  X01  „ 

oxx 

Coordinates  of  nodal  points  are  not  speci¬ 
fied. 

1XX 

Coordinates  of  nodal  points  are  specified. 

2XX 

No  element  is  associated  with  the  cor*» 
responding  nodal  point. 

Table  II  Value  and  Meaning  of  Symbolic  Word  BCCODE 


Value 


Meaning 


Adiabatic  (or  no  boundary  condition 
specified) 

Temperature  specified 
Flux  specified 
Convection  specified 


Table  HI  Value  and  Meaning  of  Symbolic  Word  IJCODE 


Value 


Meaning 


Boundary  condition  is  applied  on  segment 

(I,  J-H) 

(5) 

©  . 

(I.J-l) 


The  input  data  dsck  is  shown  in  Fi«, 
card  format  is  given  bolow. 

Card  1  TITLE  (I2A6) ' 

Coi  1-72  Any  alphanumeric  statement 

Card  2  Initial-Temperature  Card  (F10. 5) 


16,  and  the 


Col  1-10 


Initial  Temperature 


Card  3  Boundary-Temperature  Card  (4F10. 5) 
Col  1-10  Temperature  for  Side  I  =  IMIN 
11-20  Temperature  for  Side  I  =  IMAX 
21-30  Temperature  for  Side  J  =  1 
31-40  Temperature  for  Side  J  -  JMAX 

Card  4  Boundary- Flux  Card  (4F10. 5) 

Col  1-10  Flux  for  Side  I  =  IMIN 

11-20  Flux  for  Side  I  =  IMAX  l 

21-30  Flux  for  Side  J  =  1 
31-40  Flux  for  Side  J  =  JMAX 

Card  5  Boundary-Convection  Card  (8F10. 5) 

Col  1-10  Film  coefficient  >  Side 

11-20  Environment  temperature  J 
21-30  Film  coefficient 


Input  for  speci¬ 
fied  temperature 
boundary  condi¬ 
tion.  Zero 
otherwise. 


Input  for  speci¬ 
fied  boundary 
flux.'  Zero 
otherwise. 


I=IMIN 


I=IMAX 


31-40  Environment  temperature  J 
41-50  Film  coefficient 

Side 

51-60  Environment  temperature  ) 


Input  for 
convec¬ 
tive 

boundary 

condition. 

Zero 

other¬ 

wise. 


61-70  Film  coefficient 

71-80  Environment  temperature. 


J=JMAX 


I 
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TIKE  GmiSU 


FIG.  16.  DATA  BECK  FOU  AMG042 


Cards  6 
Gol 


Card  7 

Cards  8 
Col 


Card  9 
Col 

Cards  1C 
Col 


Card  11 


Element  Property  Cards  (5X,  10,  3FIQ. 0} 

I- 5  Blank 

6-I0  Identifying  number  fear  data  set*  Range,  is  from 
0  to  20.  0  is  interpreted  internally  as  1, 

II- 20  Conductivity  in  x  direction* 

21-30  Conductivity  in  y  direction* 

31-40  Product  of  density  and  specific  heat. 

Blank  Card 


Nodal-Point  Cards  (A5,  12,  13,  12,  13,  4F10.5,  315) 

lf*5  Y/ord  -  POLAR  for  polar  coordinates,  blank 
otherwise . 


6-7 
8-10 
11-12 
13-15 
16-25 
26-35 
36-45 
46-55 
56  60 
61-65 
65-70 


11 
J! 

12 
J2 

XI  (or  Rl) 
Y1  (or  01) 
X2  (or  R2) 
Y2  (or  02) 


} 

} 

} 

} 


lowest  1  (or  J)  for  line  segment. 

largest  I  (or  J)  for  line  segment. 
Zero  if  only  a  point  is  input. 

coordinates  for  lowest  I  (or  J). 
coordinates  for  largest  1  (or  J), 


TYPE  (see  Table  l). 
BCCODE  (see  Table  II). 
IJC  ODE  (see  Table  IH). 


End  Card  (15) 
1-3  END 


TimcCards  (F10.3,  F5.  0,  F5.0) 

1-10  TMAX, 

10-*15  Number  of  time  steps  from  T  to  TMAX. 
15-20  TOUT  (Prinls  temperatures  of  T  5:  TOUT). 
Blank  Card, 


b. 


Program  AMG042  outputs  in  printed  form  coor¬ 
dinates  of  each  social  point,  the  word  CODS,  for  each  nodal  point, 
the  value  of  the  boundary  condition  and  nodal  points  for  each  type  of 
boundary  condition,  element  property  data,  and  the  temperature  at 
each  nodal  point  for  each  time  requested.  Tape  Unit  6  prepares  a 
tape  which  can  be  used  as  input  to  Program  AMGQ42P. 
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Plot  Program  AMG042P 


To  assist  in  interpreting  the  output  from  program  AMG042, 
a  subsidiary  program  AMG042P  was  devised.  This  plot  program 
used  a  magnetic  tape  prepared  by  AMG042,  along  with  conventional 
card  data  input,  to  perform  the  following  functions. 

(1)  The  array  of  nodal  point  temperatures  is  scanned, 
and  the  coordinates  of  points  on  selected  temperature  contours 
are  determined  by  linear  inverse  interpolation.  The  coordinates 
are  printed  and  also  put  on  tape  for  use  with  the  plotter. 

(2)  The  values  of  temperature  along  any  specified  coor¬ 
dinate  line  are  calculated  and  printed  and  also  put  on  a  tape  for  use 
with  the  plotter. 

The  user  has  the  option  of  obtaining  all  of  the  above  information 
in  printed  form  and/or  having  these  pj  otted  by  the  Electronic 
Associates,  Inc.,  3440  Dataplotter.  The  title  is  written  at-the  bottom 
of  the  plot  sheet  (30-in.  x  30-in. )  beginning  four  letter  heights  from 
the  bottom.  Allowance  must  be  made  for  this  in  specifying  the  board 
and  data  offsets  on  the  plotter  control  card  (Card  2).-  Each  plot  re¬ 
quires  a  separate  sheet. 

There  are  some  items  on  the  plotter  control  card  (Card  2), 
such  as  board  offset  and  data  offset,  which  are  difficult  to  explain 
briefly.  Those  who  use  this  plotter  will  find  sufficient  explanation 
in  the  plotter  literature.  Those  who  have  another  plotter  available 
can  adapt  this  program  to  their  particular  needs.  If  no  plotter  is 
available,  the  tape-writing  instructions  should  be  removed  from  the 
program  and  the  printed  output  can  be  used  for  manual  plotting. 

Only  IPRINT  need  be  specified  on  Card  2  if  a  plotter  is  unavailable. 

a.  Input  Data 

The  following  input  data  must  be  included  along  with 
the  tape  created  by  logical  unit  6  in  the  temperature-calculation 
program  AMG042,  AMGQ42P  uses  logical  unit  11  to  read  the  input 
tape  and  writes  an  output  tape  on  logical  unit  12, 

The  input  data  deck  setup  is  shov/n  in  Fig.  17  and 
the  card  format  is  given  below. 


Card  l 
Col 
Card  2 
Col 


*  Card  3 
Col 


Card  4 
Col 


TITLE  (12A6) 

i~72  Any  alphameric  statement 
Plotter  Control  Card  (4F5.2,  3F10.5,  15) 


I- 5 
6-10 

II- 15 
16-20 
21-30 
31-40 
41-50 


SC  LX 

Size  factor  in  x  direction. 

SCLY 

Size  factor  in  y  direction. 

BOFFX 

Board  offset  in  x  direction. 

BOFFY 

Board  offset  in  y  direction. 

DOFFX 

Data  offset  in  x  direction. 

DOFFY 

Data  offset  in  y  direction. 

SL 

Letter  height  in  inches  (0  if  no  letter, 
ing  is  desired). 

51-55  IPRNT  Print  control.  0  =  print.  1  =  no 

print. 

Time  Control  Card  (1A5,  9E8.  l) 


1-5 

6-13 

14-21 


W1 

TIME(l)* 

TIME(2) 


"TIME=<» 

Time  values.  Only  TIME(l)  may  be 
zero.  A  large  positive  value  of 
)  time  will  skip  to  the  next  problem. 

. *  ‘  A  negative  value  of  TIME(l)  will  end 

70-77  TIME(9)  the  program  and  rewind  tapes. 

Temperature  Contour  and  Coordinate  Card  (1A5,  9E8. 1) 


1-5 

6-13 

14-21 

70-77 


W2 

TEMP(l)  or  COORDINATE(l) 
TEMP(2)  or  COORDINATE(2) 


TEMP(9)  or  COORDINATE^} 


TEMP=,  X=,  Y=, 
R=,  Z=,  or  T=. 
(TEMP=  for  tem¬ 
perature  contours. 

►  X=,  Y=,  R=,  Z=,  or 
T=  for  temperature 
on  constant  coor¬ 
dinate.)  Must  have 
a  Card  5  if  coordi- 
J  nate  is  given. 
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Card  5  Coordinate  Scale  Card  (2F  10.5) 

(Use  this  card  only  when  W2  is  X=,  Y=,  R=,  Z=,  or 
T=. ) 

Col  1-10  TNORM  Normalizing  value  for  temperatures 
(Height  of  temperature  scale  in  inches)  s 

10  X  (SIZE  T)  x  ^ 

(TNORM) 


COORDINATE  SCALE  CARD 


TEMPERATURE  CONTOUR  AND 
COORDINATE  CARD 


TIME  CONTROL  CARD 


PLOTTER  CONTROL  CARD 


TITLE  CARD 


FIG,  »7.  DATA  DECK  FOR  AMG042P 
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111.  ...TRACT 


A  new  numerical  method  for  the  solution  of  heat  conduction  problems  in 
thermally  anisotropic,  nonhomogeneous  bodies  of  complex  geometry  was 
devised  which  is  based  on  a  discretization  concept  developed  i a  the  matrix 
analysis  of  structures.  This  discretization  method,  commonly  referred  to 
as  the  finite  element  method,  reduces  the  problem  formulation  to  the  solution 
of  a  matrix  equation  for  the  nodal  point  temperatures  of  the  assembly  of 
finite  elements.  The  resulting  matrix  equation  is  stable  for  any  time  step. 
The  method  is  extremely  flexible  and  easy  to  apply.  The  method  was  applied 
by  writing  a  computer  program  for  the  solution  of  heat  conduction  problems 
in  plane,  thermally  anisotropic,  nonhomogeneous  bodies. 
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